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Abstract. In terms of statistical fluctuations, stellar population synthesis models are only asymptotically correct 
in the limit of a large number of stars, where sampling errors become asymptotically small. When dealing with 
stellar clusters, starbursts, dwarf galaxies or stellar populations within pixels, sampling errors introduce a large 
dispersion in the predicted integrated properties of these populations. We present here an approximate but 
generic statistical formalism which allows a very good estimation of the uncertainties and confidence levels in any 
integrated property, bypassing extensive Monte Carlo simulations, and including the effects of partial correlations 
between different observables. Tests of the formalism are presented and compared with proper estimates. We 
derive the minimum mass of stellar populations which is required to reach a given confidence limit for a given 
integrated property. As an example of this general formalism, which can be included in any synthesis code, we 
apply it to the case of young (t < 20 Myr) starburst populations. We show that, in general, the UV continuum 
is more reliable than other continuum bands for the comparison of models with observed data. We also show 
that clusters where more than 10 5 Mq have been transformed into stars have a relative dispersion of about 10% 
in Q(He + ) for ages smaller than 3 Myr. During the WR phase the dispersion increases to about 25% for such 
massive clusters. We further find that the most reliable observable for the determination of the WR population 
is the ratio of the luminosity of the WR bump over the H/3 luminosity. A fraction of the observed scatter in the 
integrated properties of clusters and starbursts can be accounted for by sampling fluctuations. 
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1. Introduction and motivation 

In the past few years, the increasingly detailed observa- 
tions of stellar populations in a wide variety of environ- 
ments - from stellar clusters to high redshift galaxies - 
has driven the development of increasingly complicated 
evolutionary synthesis codes. A great care has been given 
to the use of updated physical ingredients (like more re- 
alistic model atmospheres, grids of observed spectra, stel- 
lar tracks covering more evolutionary phases, increased 
wavelength range, etc) and technical details (interpolation 
methods, convolutions, etc), leading to a situation where 
the outputs produced by different codes are roughly con- 
sistent with each other, although they perhaps disagree 
on the interpretation of a given observable. Systematic 
internal errors can be addressed comparing the outputs 



resulting from different inputs or assumptions in a given 
code (e.g., Bruzual 2001), while the cross-comparison of 



different codes with the same inputs gives an estimate of 
the external errors (e.g., Lcithcrer et al. 1996). 

In contrast, the study of statistical errors has received 
comparatively little attention, even though it seems likely 
that the observed scatter in the observed properties of 
systems with a relatively small number of stars (stellar 
clusters, starbursts, dwarf galaxies or stellar populations 
in pixels), can be accounted for by sampling fluctuations 
from a given (and perhaps universal) stellar initial mass 
function (IMF). 
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In this context, Barbaro & Bertelli (1977) pioneered 
the estimation of the dispersion of observables in stel- 
lar clusters by computing Monte Carlo simulations con- 
strained to contain a given number of post- main se- 
quence (PMS) stars (so as to reproduce the number of 
observed PMS stars in the colour-magnitude diagram of 
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a given cluster). They concluded that under the assump- 
tion of a Salpeter IMF, the observed scatter in the inte- 
grated colours of young open clusters could be reproduced 



by their models. Updated analyses, by Chiosi, Bertclli & 



Bressar (1998), Girardi & Bica (1992) and Brocato et al 
(1999^ among others, yield the same results: the disper- 
sion in a colour-colour or colour-magnitude diagram can 
be accounted for by these stochastic effects, at least in 
part. 



of the most relevant observables, and apply it, as an illus- 
trative example, to the case of young stellar populations 
(ages smaller than 20 Myr), produced by statistical fluc- 
tuations in the IMF. The effect of Poissonian statistics on 
time-integrated observables such as supernovae rates, in- 
tegrated mechanical energy output and chemical yields, as 
well as major problems in the time interpolations used in 



synthesis codes, are addressed in a third paper (Cervine 



At older ages, the presence of a few cool but very lumi- 
nous stars produces the same effects. In the near-infrared, 
the small number of thermally pulsating AGB stars can 
produce strong fluctuations in th e JHK band s, as pointed 



et al. 2001b). Improvements in track interpolation tech- 



out by Santos fc Frogcl (1997 ) (see also Langon and| 



Mouh qine 2000| , for a dissenting view). Even in the op- 
tical bands, the integrated colours of globular clusters can 
be strongly influenced b y the presence of post AGB stars, 
as 



Brocato et al. (1999 ) have shown 



A similar problem arises when dealing with stellar 
populations in pixels, particularly in images or frames of 
nearby galaxies where each pixel could include a limited 
number of stars which however account for most of a given 
observable. See Buzzoni (1993) and Rcnzini (1998| ) for fur- 
ther details. 

In the context of starbursts, Mayya (1995| ) performed 
a similar analysis, this time constraining the number of 
main sequence stars with spectral type earlier than 04, 
finding that while the errors in colours were larger when 
red supergiants appeared, they decreased when the mass 
of the cluster increased, being negligible for clusters more 
massive than about 10 5 M . 

There is, however, no systematic study of the effects 
of these sampling fluctuations on the integrated proper- 
ties (such as colours or equivalent widths of emis sion and 



absorption lines) of simple stellar populations. Buzzoni 
(1989|)]pioneered this field by presenting a simple formal- 
ism based on Poissonian errors, and applied it to synthesis 
models of old stellar populations. Although powerful, this 
formalism does not account for some important sources of 
errors and is not entirely rigorous. 

While it is clear that the foreseen improvements in 
synthesis codes will attempt to refine their realism, the 
detailed comparison with observed properties will still suf- 
fer from the limited number of stars which contribute to 
some observables. The present paper is the second of a se- 
ries whose objective is to study the accuracy of evolution- 
ary synthesis models when compared with observational 
data, and to point out the oversimplifications used in such 
models. This on-going project reviews current problems 
in synthesis codes and shows possible ways to overcome 
them, providing a guideline to include improvements in 
any synthesis code. 



The first paper (Cervino et al. 2000) pointed out the 
problem of confidence limits in synthesis models set by 
the presence of statistical fluctuations in the IMF, where 
the fluctuations were estimated by Monte Carlo simula- 
tions. In this work we present a statistical formalism to 
estimate quantitatively the fluctuations expected in some 



niques will be discussed in a forthcoming paper. We will 
complete the series by extending the results presented in 
this paper to different star formation scenarios. 

The structure of the paper is the following. In Sect. 
H we present our statistical formalism, which can be im- 
plemented in any synthesis code. In Sect. [| we use the 
formalism to derive the dispersion of some relevant ob- 
servables of young star-forming regions, highlighting the 
importance of a proper use of covariance factors in the 
evaluation of the dispersion. In Sect. || we compare the 
analytical dispersion and confidence levels obtained with 
our formalism with the ones obtained via Monte Carlo 
simulations. We discuss the results and draw our conclu- 
sions in Sect. ^. In Appendix [A| we detail error propa- 
gation in observables, giving practical examples. All the 
results of the paper are available from our web server at 
tittp : //www. laef f . esa. es/users/mcs. 



2. On sampling and Poissonian fluctuations 

Let us assume that TVtot stars are observed, with masses 
distributed between m\ ow and m up . The mass rrii of the i- 
th star is a random variable whose probability distribution 
function is given by §(rrii) the stellar initial mass function. 
That is. 



$(m) dm — 1 



(1) 



For the case of a single power law IMF in that mass 
range, $(m) ~ m~ a , so that masses can be generated by a 
simple transformation from a random variable u uniformly 
distributed in the interval [0, 1] with 



m = [(1 - u)ml ow a + 



u m„ 



L _ al l/(l-a) 



(2) 



and similarly for piece-wise power laws in different mass 
ranges. Codes which use Monte Carlo sampling of the IMF 
must follow the evolution of each star generated in this 
way. 

The number JV, of stars of mass mi is another ran- 
dom variable, whose distribution function follows Poisson 
statistics. That is, the only parameter of the Poisson dis- 
tribution is precisely the value of the IMF at that mass 



dm = dNi I N tot = $(m,) dm, 



(3) 



since obviously the total number of stars is J dNi = N to t- 
Note that stellar evolution does not change the Poissonian 
nature of Nt , and that the normalised random variable n, 
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Mass [M Q ] 

Fig. 1. Ratio of variance over mean value as a function of 
mass for 1000 Monte Carlo simulations with Ntot = 10 3 
stars (upper panel) and N to t = 10 4 stars (lower panel). 
Solid line: a 2 (rii)/ni using a mass-bin of 0.5 M Q for the 
entire mass range. Triangles: a 2 (ni)/rii using a mass bin 
of 0.5 M Q in the mass range [5, 120] M Q , and a mass bin 
of 0.01 M Q when 2M Q < m < 5M Q . 



is almost a Poisson variable, since the ratio of a Poisson 
variable w ith a constant] 1 ] is not ne cessarily a Poisson vari- 
able (e.g., |Kcndall fc Stuart 1977| ). 

A simple test to check the Poissonian nature of the 
distribution of rij is simply the ratio of its variance a 2 (ni) 
to its average value <ni> as a function of mass m^. This 
ratio should be close to unity for the entire range of masses 
considered. The situation is illustrated in Fig. |l]with 1000 
Monte Carlo simulations of clusters with 10 3 and 10 stars 
with a Salpeter IMF slope (a — 2.35) in the mass range 
2 - 120 M Q . Figure [l] shows that despite the fact that 
rii is the ratio of a Poisson variable with a constant, it 
is Poisson-distributed to within 10% even in clusters with 
as few as N tot = 10 3 stars. Note that codes which use a 
binning in mass must ensure that the size of the bins is 
small, otherwise the correlations induced may produce a 
systematic effect, especially for the lower mass stars, as 
shown on Fig. [TJ. 

The fact that rii can be approximated by a Poisson 
variable makes it possible to apply a proper statistical for- 
malism, although it is slightly more subtle than one could 
naively expect. Take for instance the integrated monochro- 



1 Formally iVtot is also a random variable (the sum of random 
variables) but we assume here that is is fixed, that is, there are 
actually only N to t- 1 independent variables. For 7Vt ot ~ 10 3 or 
more this makes a tiny difference and greatly simplifies the 
formalism. 



matic luminosity C\ of the stellar population of N to t stars 
at some wavelength A at some age t. This is simply given 
by the linear combination of individual monochromatic 
luminosities t\t of stars of mass rrii at the evolutionary 
stage of age t: 



= Ntot E n i 



(4) 



For a given age t, l\i is a constant, so that C\ is 
a weighted sum of Poisson-distributed random variables, 
and would itself be a Poisson variable if the weights would 
take integer values, since the sum of Poisson variables is 
also a Poisson variable. Hence define B A as 



B A = 1 + int[min(<? Al )] 



so that 



C\{t) = Not (1 + B A ) 



iVtot 

E 

i=l 



l + int(l Ai (t)) 
1 + B, 



(5) 



(6) 



In this case, and up to a constant, C\ would be a 
proper Poisson variable, the (integer) sum of Poisson vari- 
ables, with parameter given by the sum of the individual 
Poisson parameters. However, as long as we are interested 
in dispersions and not actual distribution functions, we 
can disregard the actual distribution function of C\ as 
long as it is a linear combination of random variables. It 
is convenient to introduce at this stage the random vari- 
able 



Wi = ruNot/Mt, 



N/Mt t 



(7) 



which gives the normalised mean value of the number 
of stars of mass nij to the total mass of the cluster 
Mtot = m i assumed to be a constant^. The integrated 
monochromatic luminosity per unit stellar mass, or spe- 
cific luminosity, is 



L A = C x I M tot 



(8) 



With this notation, the expectation value of L A at some 
age t becomes 



<L A > 



(9) 



The actual integrated specific luminosity L A will fluc- 
tuate around this average value with a variance given by 



var(L A ) = <(L A - <L A >) > 

= ^Wiilt + ^IxitxjCOviw^Wj) (10) 

i i/j 

2 Again, formally M to t is also a random variable with a mean 
value, for a Salpeter IMF between 2 and 120 M Q , of <M to t> = 
5.97iVtot and a variance va.r(M to t) = 7G.28N to t, but we will 
assume here that it is a constant to simplify the formalism. 
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where the covariance is given, as usual, by 
cov(wi,Wj) = <(wi — <Wi>) (wj — <Wj>)> 



(11) 



Note that Eq. [Tcj is an exact result since La is a linear 
function of Wi . The random variables Wi and Wj are inde- 
pendent, so their covariance vanishes. Note however that 
codes which produce bins of mass and follow the aver- 
age evolution in that bin introduce a non-zero covariance 
which must be taken into account, particularly if the bin 
size is large. We have performed several tests with dif- 
ferent mass binning, and found that when the mass bin 
is smaller than 0.5 M©, the covariance terms can be ne- 
glected for M>5 M©. Furthermore, if the mass bin is made 
less or equal to 0.1 M©, the covariance terms also vanish 
in the mass range 2-5 M©. 

So far, we have shown that La (derived from C\) is 
almost a Poisson variable, so we may expect its variance 
(j 2 (L\) to be similar to its mean. This will not be ex- 
actly so, as discussed above. There are two possibilities 
here. The first one is to assume that the Poisson distri- 
bution is a good approximation to the distribution of La . 
In this case, we can derive approximate c onfidence limits 
in the standard way (e.g., Gchrcls 1986 ) and then infer 
confidence limits for derived quantities, such as magni- 
tudes, colours, equivalent widths, etc. This proves to be 
rather cumbersome. The second possibility is to disregard 
the actual distribution function of the specific integrated 
monochromatic luminosity, and deal with variances only, 
which can then be propagated for any derived quantity 
(see Appendix A). 

In this context we can introduce, following Buzzoni 
(1989)2 ^ e conce Pt °f an effective number of stars 
N* S (L\), from the relative error in the integrated 
monochromatic luminosity 



var 



(La) 1/2 _ (Ei^hWi! 



1/2 



La J2i w ihi 



V Mot Y,i W t^i 

where we have made use of the fact that var(u;i 



(12) 



var(^)/M t i 



Wi/Mt t given that it is approximately 



Poisson. Even if the distribution function is not Poisson, 
one can still define this effective number as an approxi- 
mate Poisson dispersion. 



As Buzzoni (1989) emphasized, N* s does not corre- 
spond to a real number of stars, but rather to an estima- 
tion of the relative error, in the Poisson limit. It can be 
interpreted as a number of contributors to a given wave- 
length at that age, so that if this number is small, we may 
expect large fluctuations in the mean value of the lumi- 
nosity at this wavelength. In the case where all the stars 
have the same luminosity, £\i — l\, and since J2 w i = 
JVtot /Mtot then A* ff (LA) = N to t, and all stars contribute 
the same, obviously. In the general case, however, A* ff (LA) 
can never exceed N totl by definition. Hence small values 



indicate that the average value must be interpreted with 
caution, the variance around the mean value being as large 
as the mean value, typically. The ratio N* s (L\)/N tot can 
also be interpreted as a measure of the statistical entropy 
of the l\i values (Buzzoni 1993): for a given M tot , a larger 
value of iV* ff (LA) means a lower dispersion in 4t- We can 
also rescale N* s to the total mass in stars, following the 
convention used in population synthesis : 



A eff = N; s / M tot 



(13) 



and the definition of N* s can be generalized to almost all 
observables of a synthetic population (ratios, equivalent 
widths, colours etc., see Appendix A for more details). 
This formulation has two main advantages: 

— For a cluster of a given mass, N e s gives a simple and 
easy to estimate measure of the expected dispersion 
of a given observable. This dispersion depends on the 
number of stars in the cluster which effectively con- 
tributes to that observable. It must always be com- 
puted if a meaningful comparison with observed quan- 
tities of clusters of that mass is made. 

— The comparison of the different values of N e g for dif- 
ferent observables (say, luminosities in two different 
passbands) gives a direct estimation of the difference in 
the statistical entropy of the observables, which helps 
to establish which observable is more reliable and thus 
better constrained by observations. 

The N c g formalism can, and should, be included in any 
synthesis code, and the corresponding values of iV c ff(x) 
must be given along with every computed observable x so 
that the dispersion in x can be assessed. Before addressing 
the comparison of the dispersions computed with this an- 
alytical formalism with the ones obtained more precisely 
via Monte Carlo simulations in Sect. ||, it is worth giving 
a practical example of its application. 



3. A worked example 
starbursts 



the case of young 



The first application of a simplified version of this formal- 
ism was made by Buzzoni (1989, |1995 ) in his study of old 
stellar populations, where the presence -or otherwise- of 
some evolutionary phases affected important observables 



in the optical (see also Langon and Mouhcine (2000) for 
the case of infrared colours). Rapid evolutionary phases 
also occur in young starburst populations, and so as an 
illustration of the usefulness of the formalism, we ap ply it 
here to an updated version of the code presented in 



Mas- 



Hesse and Kunth (199l|) ; |Ccrvino and Mas-Hesse (1994J) ; 
Cervino et al. (2001a). The main updates to the code in- 



clude the following : 

1. Inclusion of the full set of Geneva evolutionary tracks 
i ncluding standard and enhanced mass-loss rates 
flMeynet et al. 1994| ). 
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2. Inclusion of metallicity-dependent atmosphere models 3.1. Continuum properties 



for normal stars from Kurucz and CoStar ( Schaerer & 



De Koter 1997) and the atmosphere models for WR 
stars from ftchmutz ct al. (1992j ). 



3. Inclusion of an analytical IMF formulation using a 
modified dynamical mass-bin from the Dec. 2000 re- 
lease of Starburts99 ( Leitherer et al. 1999| ). We have 
also kept the original Monte Carlo formulation. The 
dynamical mass-bin subroutine used here makes use of 
a more restrictive mass binning than the one used in 
Starburts99 with a higher numerical precision in order 
to ensure a narrow enough mass bin for the computa- 
tion of the dispersion. The typical size of the dynamic 
mass bin used in our code is 10 mass points. This nu- 
merical resolution does not produce however any sig- 
nificant changes in the mean properties obtained by 
both codes. 

4. Use of parabolic interpolations in time for the com- 
putation of the isochrones (see Cervino et al. 2001b, 
Paper in of this series, for a full study and justifica- 
tion). 

The code also includes an intermediate track that takes 
into account the discontinuity in the evolution of WR and 
non-WR stars. Such an intermediate track was introduced 



by Cervino and Mas-Hesse (1994 ) and will be discussed in 
Paper iv of this series. 

The complete predictions of the updated set of evo- 
lutionary synthesis models will be presented elsewhere. 
For the present example we will focus on solar metallicity 
s imple stellar popula tions with standard mass-loss rates 
(fcchallcr et al. 1992| ), a Salpeter IMF slope m a mass 
range 2 - 120 M Q , an Instantaneous Burst star-formation 
law and ages from to 20 Myr. We will use throughout 
this work an analytical approximation of the IMF, we will 
not use the intermediate tracks mentioned above, and we 
will neglect the presence of binary systems, so that our 
dispersion estimates can be compared to those obtained 
by other synthesis codes with similar inputs. 

In the following subsections we discuss the dispersion 
estimated for a few observables predicted by the code: 
continuum luminosities, V-K colour, H/3 equivalent width 
and some WR ratios. We focus here on the evolution of 
the dispersion of these observables, but we also show their 
general evolution in order to understand their dispersion. 
A detailed study of the evolution of these observables 



has been discussed elsewhere (see Cervino and Mas-Hesse 



1994) for a more detailed descr iption of the evolution 



and Leitherer fc Heckman (1995 ), for a comparison with 
other synthesis models, showing that both codes give quite 
similar results)^. 



3 The only relevant difference between the so- 
lar metallicity models presented here and the 
ones from the Dec. 2000 release of Starburts99 



( ittp : //www. stsci . edu/science/starburts99/ ) is the 

use of different atmosphere models. Incidentally, the re- 
sults presented here can be used to evaluate dispersions 



Equation |l2| can be applied in a straightforward way to 
the continuum predicted by the synthesis code. We only 
need to take into account that each l\ i is the luminosity 
of a star of initial mass mj with a relative weight given by 

Wi. 

A small complication arises in the case of young star- 
forming regions in that each massive star also contributes 
to the ionisation of the cloud around the cluster, and thus 
to the nebular continuum. This contribution depends on 
the properties of the gas in the nebula (the electronic tem- 
perature T c , the electronic density n c , and the metallicity, 
in the simplest case). We have used in the following the 
same inputs as Mas- Hesse and Kunth (1991 ): T c = 10 4 K, 
n c = 100 cm" 3 and N(He n)/N(H n) = 0.1. We have also 
assumed that only a fraction / =0.7 of the ionising pho- 
tons is absorbed by the gas. 

In this case the luminosity £ \ i at a given wavelength 
has two components, i.e. l\i = + l n> i where is the 
stellar luminosity, and is the contribution of the star 
to the nebular continuum. The scaled variance in the total 
(stellar plus nebular) specific luminosity is then 



<„ t = vartLtoOM^ 1 = Y J W i l\ i = 

i 

= ^ w i l l,i + X! W ^n,i + 2 X! W i l *,i l n,i ( 14 ) 
i i i 

The last term of the equation is the covariance between 
stellar and nebular continua for a given mass or star i : 



cov(L*, L„) = 2J Wi h,i l n ,: 



(15) 



Note also that the scaled dispersion of the specific 
monochromatic luminosity er(Ltot) has units of erg s^ 1 
A-!M- 1/2 . 

The overall evolution of the continuum from 100 A to 
10 4 A is illustrated in the upper left panel of Fig. |[ At 
the beginning of the burst, massive stars are the dominant 
source of luminosity from the UV to the near infrared. The 
extreme-UV continuum (below 500 A) appears with the 
first WR stars (at about 2.5 Myr) and lasts until the end 
of the WR phase (at about 5 Myr). The lower left panel 
shows the evolution of N e g(L\) as a function of age and 
wavelength. The figure shows that the extreme-UV con- 
tinuum produced by WR stars has a very small N e g(L\) 
value, reflecting small number statistics and the statistical 
entropy of the luminosities of these stars in the synthetic 
cluster at these wavelengths. It also shows a large disper- 
sion in the optical and near infrared fluxes at ages between 
2.5 Myr and 8.5 Myr. 

The evolution of the continuum for some selected ages 
and the corresponding N e g(Li\) values are shown on the 



in the Starburst99 predictions, when using the same input 
parameters. 
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2000 4000 ^000 SODC ' D 5000 10^ 

wavel&nqth [Al wavdenqti [41 



Fig. 2. Specific continuum flux. Top left: Evolution of the overall spectrum. Units are L Q A 1 Mg 1 . Bottom left: 
Effective number of stars for the specific flux, in units of Mg 1 , for the overall spectrum. The right panels show the 
same quantities for selected ages. 



right panels of Fig. p|. The compa rison of the lower right 
panel with Fig. 2 of Buzzoni (1989) is interesting. Buzzoni 
presents N e g(L\) for a 15 Gyr synthetic cluster, showing 
an abrupt change in N c g(h\) at the Balmer jump. In spite 
of the differences in the synthesis codes, ages and metallic- 
ity, the same pattern is found in our Fig. ^. The compar- 
ison of both figures reveals however the different nature 
of the Balmer discontinuity for stellar clusters with differ- 
ent effective temperatures. In clusters dominated by hot 
massive stars, the blue continuum is dominated by the 
nebular emission and is more luminous than the red con- 
tinuum, while clusters with cooler stars show the opposite 
behaviour. This means that the discontinuity in jV e ff (L^) 
at the Balmer jump reflects the increase in the statis- 
tical entropy of the luminosities of the stars (including 
the stellar and the nebular contributions) at this wave- 
length. At small ages (0.1 Myr), the cluster is dominated 
by hot stars, emitting shortwards of the Balmer contin- 
uum: the statistical entropy for A < A(Bac) is lower than 



that found at A > A(Bac). At 3.1 Myr, the nebular con- 
tribution to the continuum decreases rapidly. The stellar 
continuum contributes most to the total flux so that the 
statistical entropy for A < A(Bac) becomes now larger 
than at A > A(Bac). At 11.1 Myr, the remaining stars are 
so cool that the nebular emission no longer plays a role in 
the Balmer discontinuity. At older ages (cooler stars) the 
Balmer discontinuity appears completely inverted. Let us 
illustrate this in more detail by studying the correlation 
of the stellar and nebular fluxes : 

— For a 0.1 Myr cluster, the nebular emission contributes 
more than 30% to the total (stellar + nebular) flux at 
wavelengths longwards ~ 4000 A. The emission at such 
wavelengths is therefore correlated with the ionising 
flux of the stars which dominate the emission at shorter 
wavelengths, and so N c g(L\) is essentially flat. 

— At 3.1 Myr, the nebular emission contribution at 4000 
A drops to 5%. The correlation with the ionising flux 
of massive stars disappears and this is reflected by a 
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decrease in the value of N c s(L\). Also, the emission 
of the different massive star populations (i.e. the £\i 
values) is very variable and depends strongly on the 
amount of stars that become WRs. Such a situation 
produces a large dispersion in the ionising, optical and 
infrared continua. 

— The dispersion in the ionising continuum vanishes at 
ages older than 5 Myr (when the WR population dis- 
appears) , but remains large in the optical and infrared 
continua due to the presence of red supergiant (RSG) 
stars. Note that WRs and RSGs have a rapid evolu- 
tion, they are very luminous, and a small amount of 
these stars may dominate the overall continuum. As 
long as their total number is small, the intrinsic dis- 
persion at such phases of the evolution of the cluster 
will be large. 

— Finally, a 11.5 Myr old cluster shows a small Balmer 
jump, and all the stars in the cluster have "similar" 
UV spectra, i.e. N c tf(L\) reaches a maximum value at 
this wavelength. The spectrum of the individual stars 
becomes more and more different towards the red, as 
shown by the decay of N e g(L\) at longer wavelengths. 

In general, the UV continuum has a larger iVcff(L,\) 
value than the optical continuum (and thus a larger statis- 
tical entropy) . This also reflects the fact that the UV con- 
tinuum is more homogeneous along the spectral types, i.e. 
small differences in effective temperature lead to similar 
UV spectra. The optical continuum depends more on the 
effective temperature, and so small variations in the popu- 
lation produce larger effects. Such a dependence becomes 
more and more important towards longer wavelengths. 

At the other extreme of the wavelength range, in the 
ionising continuum, the situation is even more dramatic 
than at IR wavelengths. The differences in the ionising 
continua of different stars are generally quite large, and 
so this part of the continuum will be the most affected by 
statistical fluctuations, as shown below. 

3.2. Ionising continuum 

We have calculated the predicted ionising flux shortwards 
of the He+, He and H° edges, Q(He+), Q(He°), and 
Q(H°) respectively. The corresponding values as well as 
-Weff(Q) and cr(Q) are shown in Fig. [| As before, the val- 
ues shown in Figure || are normalised to the total mass 
transformed into stars M tot . 

An important quantity in this context is the mini- 
mum mass of the cluster required to ensure a dispersion 
smaller than 10% in a given observable, M 10 %. Clusters 
with smaller masses will have a dispersion larger than 10% 
in that observable. For instance, the right axis in the bot- 
tom left plot of Fig. | gives the mass of the ensemble of 
stars required to ensure a dispersion in the Q values lower 
than 10%, denoted as M 10 %(Q). Note that the relation 
between the numbers of stars present iVtot and the total 
mass of the cluster depends on the IMF and the mass lim- 




5 10 15 



Age [Myr] 

Fig. 3. Ionising photons. Top left: Evolution of the specific 
number of emitted ionising photons able to ionize He + , 
Q(He+), He , Q(He°), and H°, Q(H°), in units of photons 
s _1 Mg . Bottom left: Effective number, in units of Mg , 
for the different ionising photons (see text for definition). 
The right axis shows the minimum mass required to ensure 
a Poissonian dispersion smaller than 10% assuming a mass 
range 2 - 120 M Q with a Salpeter IMF slope. Top right: 
cr(Q) in photons s -1 Mg 0,5 units for the same quantities. 

its adopted, and so the values of M 10 % must be rescaled 
if different assumptions are made regarding these values. 

Note that the scaled dispersion cr(Q) is larger than the 
normalised value of Q. It does not mean that we can obtain 
negative values of Q, because the probability density dis- 
tribution is quasi-Poisson, i.e. the mean value and the cor- 
responding error is not Q±cr(Q). Take for instance an ex- 
treme example, a cluster with a total mass Mtot = 120Mq 
made up by stars with masses in the range 2-120 Mg at an 
age of 0.1 Myr. The predicted value of the scaled (specific) 
average of Q(H°) is 4xl0 47 photons s _1 Mq , i.e. the ac- 
tual value of Q(H°) for such a cluster would be 4.8xl0 49 
photons s _1 . But the actual Q(H°) may vary between 
photons s _1 (if none of the stars in the cluster produce the 
ionising flux) and ~ 5 x 10 50 photons s _1 (if the cluster 
is formed by only one star of 120 M Q ). In fact, the model 
predictions at the 90% confidence level for such a cluster 
is Q(H°) = 4.8 x 10 49 photons s" 1 (but see Sect. | 
for details and cautions on such an estimation). 

This figure also shows another important feature: the 
relative dispersion (the inverse of N e g) becomes larger for 
more energetic edges. It is often assumed that the uncer- 
tainty on Q(H°) is representative of the uncertainty of 
the ionising continuum at shorter wavelengths. Figure |^ 
shows on the contrary that this assumption is false, since 
the shorter the wavelength, the lower the related N e g(Q), 
and the larger the relative dispersion. 

Note that to ensure a relative dispersion of 10% in the 
flux below the Q(He + ) break more than 5 10 4 M Q must 
be transformed into stars (in the mass range 2 - 120 M© 
with a Salpeter IMF). The situation becomes even worse 
when WR stars appear in the cluster. In this case a mass 
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Fig. 4. Luminosities. Top left: Evolution of the specific 
monochromatic luminosity at 4861, 5500 and 22000 A in 
units of erg s _1 A -1 Mq 1 . Bottom left: Effective number, 
in units of Mq 1 , for the different monochromatic luminosi- 

in units of 



ties 
erg s 



Right axis as in Fig. |3|. Top right: o~(L\) 
1 A -1 M B ' 5 for the same quantities. 



of about 10 6 M Q is needed to get an uncertainty of about 
10%. Note that the corresponding N c g(Q) values are lower 
than the UV continuum ones, so that the dispersion in the 
emission line spectrum is expected to be much larger than 
in the continuum. 

In the past few years the continuum of evolutionary 
synthesis models has been extensively used in the com- 
putation of grids of photoionisation models. Although the 
general predictions of photoionisation models can be valid, 
the fitting of observed data with such grids, or the com- 
parison of the grids with a set of observed sources, must 
take into account the possible dispersion due to these sam- 
pling fluctuations in the stellar population. The expected 
dispersion on the emission line spectrum will be evaluated 
in a forthcoming paper. 

3.3. Monochromatic luminosities and colour indexes 

Figure |i| shows the evolution of the continuum luminosities 
at the wavelength of the H/3 line, at 5500 A (V band) and 
at 2.2 pm (K band). The figure also shows N e fi(L\) for 
the given wavelengths and the required minimal cluster 
masses to ensure a relative dispersion smaller than 10%, 
M 10 %(La), as well as the associated dispersion o-(L\). 

As it was pointed out above, at the beginning of the 
burst the nebular contribution dominates, and so N c f[(L\) 
is similar at all wavelengths. Later on, N c g(L\) decreases 
when the fast evolving stars (WR and RSG) appear in the 
cluster. The dispersion decreases for optical wavelengths 
when WR stars disappear. 

The re sults from this figure can b e compared with 
Tabic I in Langon and Mouhcinc (2000). Both result s are 
quite in agreement, as well as with Buzzoni (1993 ): IR 
luminosities require in general a larger amount of mass 
transformed into stars to ensure a given relative disper- 
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Fig. 5. V-K colour and its covariance. Top panel: 
Evolution of the V-K colour index for a cluster where 
the nebular contribution is included (left) and for a clus- 
ter with the stellar contribution only (right). Bottom pan- 
els: cr(V-K) in units of mag Mq 5 obtained with different 
assumptions about the covariance terms for the two sce- 
narios mentioned above. 

sion. Equivalcntly, IR luminosities are much more affected 
by the discreteness of the stellar population and are more 
prone to sample fluctuations. 

The minimal cluster masses required to ensure a dis- 
persion smaller than 10 % for the 5500 A and 2.2 pm lu- 
minosities obtained by Langon and Mouhcinc (2000] ) are 
however a bit larger than the ones obtained by us. This 
can be understood taking into account that they used a 
different mass range (0.1 - 120 M Q ) and that they did not 
include the nebular contribution in their calculations. 

In Fig. H we show the evolution of the V-K colour 
index and its (scaled) dispersion. Since the results depend 
strongly on the proper evaluation of the covariance factors, 
we have also considered the resulting V-K colour index 
for the case where the stellar contribution is the only one 
used. 

Note that the normalisation units of cr(V-K) are mag 
Mq 5 rather than mag Mq ' 5 . As shown in Appendix A, 
this is the case for the dispersion in ratios : the larger the 
mass, the lower the dispersion. The normalised ct(V-K) 
value at 20 Myr is about 40 mag Mq 5 , so that an intrinsic 
dispersion of « 0.4 mag is typical for a 10 4 Mg cluster. 
It increases to w 1.26 mag for a 10 3 Mg cluster (see also 
Sect J). 

An important factor in the calculation of the disper- 
sion is the covariance term. Consider the definition of the 
(linear) correlation coefficient of two quantities p(x, y) as 



p{x,y) 



cov(x, y) 



(16) 



It varies between —1 and +1 where the sign indicates 
the sense of the correlation, so that when p(x, y) is equal to 
1 the quantities are completely correlated and if there 
is no correlation. For one given star, as pointed out by 
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Fig. 6. Correlation coefficients. Evolution of the correla- 
tions coefficients for U and B, B and V, and V and K for 
a cluster where the nebular contribution is included (left) 
and for a cluster with the stellar contribution only (right). 



Buzzoni (1993 ), the luminosities in different bands are 
completely correlated since they are produced by the same 
star. But the situation changes when several stars, instead 
of an individual one, are considered, i.e. the case of popu- 
lation synthesis. This is shown in the lower panels of Fig. 

where the corresponding tr(V-K) values are shown for 
different assumptions made about the covariance. The as- 
sumptions made are (1) no inclusion of covariance terms 
(or equivalently cov(V,K) — and so p(Y, K) = 0), (2) a 
complete correlation, p(V,K) = 1, and (3) the case where 
the covariance is obtained from the model adding properly 
the contribution of each population of synthetic stars. The 
figure clearly shows that the use of p(V, K) = leads to 
an overestimate of the error while the use of p(V, K) = 1 
produces an underestimate. Note that this trend is inde- 
pendent of the presence or otherwise of the nebular con- 
tinuum. 

Lancon and Mouhcine (2000) argue that a complete 
correlation in stellar clusters can only be guaranteed for 
nearby continuum bands. We have tested this assumption 
in Fig. |^, where the correlation coefficient of some colour 
indexes have been plotted. 

The figure shows that, as expected, p(V — K) deviates 
from unity, which means that the covariance term must 
always be taken into account. It also shows that a similar 
situation happens with p(U — B) and p(B — V). The ef- 
fects become even more severe where the nebular emission 
is included (especially for U-B, since the U filter contains 
the Balmer jump). Therefore the assumption of complete 
correlation (p(x,y) = 1) produces an underestimate of 
the dispersion even in close continuum bands. This result 
highlights the importance of the proper inclusion of the 
covariance term at all wavelengths. 

3.4. Hf3 equivalent width and WR ratios 

Figure ^ shows the evolution of the EW(H/3), and the 
associated iV" cff (EW(H/3)) and <r(EW(H/3)). As in the case 
of colours, the scaled cr(EW(H/3)) is obtained dividing the 
square root of the variance by the square root of the total 
mass transformed into stars. 

The relative dispersion in EW(H/3) is similar to the 
dispersion in the continuum at such wavelengths (c.f. Fig. 
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Fig. 7. Equivalent width of H/3. Top left: Evolution 
of EW(H/3) Bottom left: Effective number, in units of 
Mq 1 , for EW(H/3). Right axis as in Fig. | Top right: 
cr(EW(H/3)) in units of A M^ 5 . 

[|) except for the very early phases of the burst. From to 
2 Myr L(486lA) and Q(H+) are highly correlated via the 
influence of the nebular continuum and it gives raise to a 
lower dispersion on EW(H/3). 

Finally the left panel of Fig. || shows the evolution of 
the blue WR bump at 4686A and H/3 luminosity ratio, 
the WR over (WR+O) and WC over WR ratios, as well 
as the corresponding N B s and a values. The right panel 
shows similar quantities for the equivalent widths of the 
4686 A WR band and 5808 A lines. 

The figure shows that the L(WR bump) /L (H/3) ratio 
is more reliable than the equivalent widths of the WR 
band and 5808 A line to determine the presence of a WR 
population in starburst galaxies (where WR population 
ratios can not be obtained directly). This is due to the fact 
that the WR population itself produces ionising photons, 
and so it is correlated to the H/3 luminosity as well. 

4. Comparison with Monte Carlo simulations 

The examples in the previous section show the usefulness 
of the formalism at evaluating, in a simple way, the dis- 
persion in integrated properties, and the insight that they 
provide on the effects of sampling fluctuations in observ- 
ables. The question that remains to be checked is whether 
the estimates of the dispersion provided by the formalism 
are realist. 

This issue can be solved performing Monte Carlo sim- 
ulations where the actual distribution function of a given 
observable is computed and its dispersion evaluated. As 
an example we show here simulations for the V-K colour 
index and EW(H/3). This allow us to test to which extent 
the formalism and the influence of the covariance terms 
are correct. 

The simulations have been done using 1000 clusters 
with A^tot^lO 3 stars, 500 clusters with 10 4 stars and 100 
clusters with 10 5 stars. In each set we have obtained 
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Fig. 8. Left Panel: WR ratios Top left: Evolution of the WR bump over H/3 luminosities ratio, WR over (WR+O) 
ratio and WC over WR ratio. Bottom left: Effective number, in units of Mq 1 , for the given observables. Right axis as 
in Fig. ||. Top right: a in units of Mq 5 (Note the normalisation units). Right Panel: WR equivalent widths Top 
left: Evolution of the EW(WR) for the 4686A band and 5808 A line. Bottom left: Effective number, in units of M" 1 , 
for the given observables. Top right: er(WR bump) and <r(5808 A line) in units of Mq 5 . 



the dispersion cr* lus (V-K) and <r* lus (EW(H/3)), as well as 
JV e * ff (EW(H/?)). The cr* lus values have been divided by 
A^ 5 , to obtain a normalised o" c ius value for each set. 
A* ff (EW(H/3)) has been divided by A to t to obtain a scaled 
JV c ff(EW(H/3)) value which can be compared with the an- 
alytical results^. 

The top left panel of Fig. || shows the comparison of the 
mean V-K colour index for the different simulations and 
the analytical solution. The values are clearly consistent 
with each other, and show not only that the analytical 
formalism is a very good approximation, but also that the 
formalism can be used in place of Monte Carlo simulations. 

There are, however, some deviations in the 10 3 stars 



cluster simulations. As shown in Cervino and Mas-Hesse 
(19940 the V-K colours at the ages where the major dis- 
crepancies appear are mostly dominated by the presence 
of RSGs. The differences between the analytical estima- 
tion of the dispersion and the real one (as obtained from 
Monte Carlo runs)is perhaps due to the low number of 
stars responsible of the colours. For such a small num- 
ber of stars the sampling is too poor to be meaningful. In 
these low mass clusters the analytical dispersion provides 
only a first order approximation to the real one, that can 
only be obtained from Monte Carlo simulations. The lower 
left panel shows the comparison of the normalised values 
of <T c i us (V-K) with the analytical value. It shows that the 
analytical dispersion and the Monte Carlo results are very 
similar. 

We also show in Fig. || the 90% confidence levels of 
the different simulations. It is important to highlight that 



4 Note that in the simulations a fixed total mass for the clus- 
ter cannot be achieved, since it is a random variable. Hence, a 
renormalisation to N c h= 1 star is used rather than to M to t = 
1 Mq (i.e. we have used a factor 0.17 = iVtot /<M to t> for the 
conversion) . 



these confidence level bands are not symmetric with re- 
spect to the mean value. 

It is also possible to obtain confidence level bands an- 
alytically, as shown in Fig. [l(]. The analytical confidence 
levels have been obtained using the approximations found 
by Gehrels (1986) for Poisson statistics. We have used 
the corresponding denormalised N c g, (i.e. N* s ) for the 
EW(H/3) and the Lv/Lk observables as the Poisson pa- 
rameter. The resulting confidence levels have been mul- 
tiplied by A/N* s with A corresponding to EW(H/3) or 
Lv/Lk- In the case of V-K, the resulting confidence level 
of Lv/Lk has been transformed into the corresponding 
V-K value. 

The analytical estimation of confidence levels for clus- 
ters with masses larger than 10 3 stars provides a very 
good approximation to the one obtained in the Monte 
Carlo simulations, showing, a posteriori, that the actual 
probability distribution function is well approximated by 
a Poisson distribution. 



5. Discussion and conclusions 

In this paper we have reviewed the effects of sampling 
fluctuations in evolutionary synthesis codes. We have pre- 
sented an approximate, and yet accurate, analytical for- 
malism which allows the estimation of Poisson dispersions 
for any observable synthetised by population codes. It is 
important to take into account that these dispersions vary 
from observable to observable and cannot be extrapolated. 
We provide some simple cases (Appendix A) for a variety 
of observables. The comparison of these analytical disper- 
sions with proper Monte Carlo simulations shows that the 
formalism is valid, and hence bypasses the need to per- 
form Monte Carlo simulations to evaluate the dispersion 
in any given observable. The formalism is simple enough 
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Fig. 9. Comparison with Monte Carlo simulations: Left Panel: V— K color index Top left: Evolution of the V-K 
color index for the analytical and Monte Carlo simulations. Bottom left: Comparison of the normalised value of er(V-K) 
for the analytical and Monte Carlo simulations. Middle and Right Panels: EW(H/3). Top middle: Evolution of 
EW(H/3) for the analytical and Monte Carlo simulations. Bottom middle: Normalised N e g (EW(H/3)) for the analytical 
and Monte Carlo simulations. Top Right: Comparison of the normalised value of cr(EW(H/3)) for the analytical and 
Monte Carlo simulations. 
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Fig. 10. Comparison of analytical confidence levels with Monte Carlo simulations. Left Panel: V-K colour index. Right 
Panel: EW(H/3) 



that it can (and must) be implemented in all the synthesis 
codes to assess the predicted dispersions in the computed 
observables. We also stress the importance of evaluating 
the covariance terms when estimating the Poissonian dis- 
persion. 

As an illustrative example of the formalism, we have 
implemented it for the case of young starbursting clusters. 
We find that in general the UV observables are the best 



defined ones (i.e. they have the lowest Poissonian disper- 
sions). The less reliable continuum range corresponds to 
the hard ionising flux. In particular, a relative dispersion 
smaller than 10% for the hard ionising flux can only be 
achieved for clusters more massive than 5 10 4 M . The 
situation becomes even worse when WR stars appear in 
the cluster. In that case a mass of about 10 6 M Q is needed 
to ensure an uncertainty of about 10%. The implications 
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of such effects on the use of evolutionary synthesis models 
as input for photoionisation ones remain to be evaluated. 

In the case of WR stars, we have shown that the most 
reliable observable for the determination of the WR pop- 
ulation is the ratio between the WR bump at 4686 A and 
the H/3 luminosities. We have also shown that the pres- 
ence of WR stars is a quite poorly-defined property of 
star forming regions, due to small number statistics. In 
particular, it is possible that the general failure of synthe- 
sis models to predict the WR population in our Galaxy 
could be solved taking into account the dispersion in real 
clusters due to Poissonian fluctuations in addition to the 
dispersion introduced by the star formation history. 

We have also shown that the dispersions and confi- 
dence levels derived analytically are in good agreement 
with the dispersions obtained from Monte Carlo simula- 
tions for clusters containing more than 10 3 stars, i.e. more 
than 5 x 10 3 M© transformed into stars in the mass range 2 
- 120 Mq following a Salpeter IMF. The confidence levels 
and dispersions obtained can be used as a first order evalu- 
ation of the number of Monte Carlo simulations needed to 
obtain reliable results, or used directly in the comparison 
with observed data. 

Acknowledgements. We acknowledge the pioneering work by 
Alberto Buzzoni which was the inspiration for this study. This 
research has made use of the LEVEL5 utility at NASA/IPAC 
Extragalactic Database (NED) which is operated by the Jet 
Propulsion Laboratory, California Institute of Technology, 
under contract with the National Aeronautics and Space 
Administration. MC has been supported by an ESA post- 
doctoral fellowship. JMMH has been partially supported by 
Spanish CICYT grant ESP95-0389-C02-02. 

Appendix A: Observables, error propagation and 
the effective number of stars iV e fr 

To illustrate the use of the N e g formalism for a variety of 
astrophysical observables, we give here a few examples to 
help the reader. For a complete review on error propaga- 
tion, see Kendall fc Stuart (1977 ) or a simplified version 



by Leo (1992) and its adaptation on the NE D server at 
tittp : //nedwww . ipac . caltech. edu/level5. 

Let us assume a random variable u which is a function 
of two random variables u — f(x,y). Let us also assume 
that ai. and c\, are the variances of these variables. The 



variance of the random variable 
mated as 



o~\, can be approxi- 



The correlation coefficient varies between —1 and +1 
where the sign indicates the sense of the correlation. 

We can give a more general definitio n of the effectiv e 
number of stars 7V e ff, as introduced by Buzzoni (1989 ), 
to any observable. As it was pointed out in Sect. 2, we 
can define a iV c ff (A) of any synthesized quantity A which 
is the sum of the contributions from individual stars (or 
populations), A = ^uJiOi, and thus it has a variance 



var(A) = ^var(wi)af (c.f. Eq. A. 4). The relative error 
is then 



var^) 1 ^ _ (X>ia?) 



1/2 



1 



(A.3) 



Note that if A is normalised to the total mass, N c g(A) 
is also normalised. So N c g(A) gives an estimate of the 
relative error (and hence, of the absolute error) for any 
possible value of the mass transformed into stars. 

A.l. Error in the multiplication by a constant 

Let us assume an observable random variable u(x) = ax 
where a is a constant, with Eq. A.l one obtains 



2 2 2 

<r„ = a u, r 



and the relative error 
1 



1 



u x ^N cS {x) y/N cS {u) 



(A.4) 



(A.5) 



This is the more general case used in synthesis codes, 
as it is showed in Eq. |lj. 

A.2. Error of a sum 

Let us now assume an observable random variable u = 
x + y. In this case 



°l = °l + tf + 2 cov(x, y) 



(A.6) 



The value of cov(x, y) depends on the observables we 
use, as shown by Eq. [l4|. 

An interesting case comes from the computation of the 
luminosity in a band, for example Lb, obtained by the 
integration of the monochromatic flux, /(A), over the filter 
response of the band, R(X): 



df 
dx 



df 
dy 



2cov(x, y) 



d A d A (A.l) LB=y'/(A)*i?.(A)dA~^/(A l )i?A,AA l (A.7) 



dx dy 



where cov(x, y) is the covariance. This equation is exact 
if the function / is linear in the x and y variables. The 
covariance can be related to the linear correlation coeffi- 
cient, p(x,y), defined as 



p(x,y) 



cov(x,y) 



(A.2) 



where we have approximated the value of the integral by 
the resulting value from the use of the trapezium rule 
and the individual monochromatic fluxes, f(Xf) and band 
transmission coefficients, R\ i . The corresponding variance 
due only to the error in the monochro mati c luminosities, 
o~f(\A = o~i, can be obtained from Eq. A.l and using Eq. 
] and Eq. O: 
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[Rxt AA, 



+ V 2 cov (/(A,), /(A,)) i? Ai i? Aj AAjAAj (A. 



i<j 



Note also that the resulting variance (and correspond- 
ing iV e ff ) must be similar to the one obtained if the syn- 
thesis code computes Lb, summing up the individual Lb ; 
and using Eq. A. 3. 



N c g{u) 
1 



N eS (x) N cS {y) 



cav(x,y) 
xy 



(A.15) 



This is the case, for example, of the computation of the 
dispersion in the mechanical energy or mechanical power 
as given by the mass-loss of the stars in the cluster and 
their terminal velocities. 



A. 3. Error of a difference 

Let us assume an observable random variable u = x — y. 
The corresponding variance is 



°l + °l - 2cov(ir, y) 



(A.9) 



A typical example in synthesis codes is the computa- 
tion of colour indexes. As an example, for the U-B colour 
index one obtains (see also below for the error of a loga- 
rithm) : 

^u-B = a l + °i - 2cov(U, B) (A. 10) 

This leads to <?u_ B — (cu — ub ) 2 if the correl ation 
coefficient is 1, as it is assu med by Buzzoni (1993] ) and 
Langon and Mouhcinc (2000 ) for colours where the corre- 
sponding fluxes originate from the same population (but 
it is not the general case for integrated colours from dif- 
ferent stellar populations). Note that the error in colour 
indices may be lower than the error in the corresponding 
colours, depending on the covariance term. 

A. 4. Error of a product 

Let us assume an observable random variable u = xy. In 
this case 



°i = x 2 a 2 + y 2(J l + 2xy cov(x, y) 
so that 

°l _ °l , G l , o cov(a;,y) 



9 2 2 

u z x y 



xy 



(A.ll) 



(A.12) 



We can rewrite the relative error as a function of N e s- 
When the correlation coefficient is zero (cov(x,y)=0) : 



1 



1 



1 



u 2 N eS (u) N cS (x) N cS (y) 



(A.13) 



while if the correlation coefficient is 1, cov(x,y)=a x a y 

(A.14) 



and 
2 1 



1 1 

+ 



In the general case: 



A. 5. Error in a ratio 

Let us assume an observable random variable u = x/y 
where y is non zero. The resulting equation for the relative 
error is 



in. ~ ln\ + _v _ 2 

u 2 x 2 y 2 



cov(x,y) 



xy 



(A.16) 



which is similar to Eq. A.12 but changing the sign which 
corresponds to the covariance term. 

In this case the normalisation is very important. Since 
the normalised quantity is cr 2 , then x 2 /a 2 , y 2 /(7y and 
xy/cov(x, y) have the same normalisation as x or y, and so 
is u 2 /a^ which has the same normalisation as x or y. And 
since u is dimensionless l/<r 2 is the normalised quantity 
(and not a 2 as in the general case). Note that the larger 
the mass of the cluster the larger should the absolute error 
be, which makes no sense. It also shows the advantage of 
using N c g instead of a 2 for the description of errors: N e s 
always keeps the same normalisation. 

Ratios are used a lot in synthesis codes (equivalent 
widths, populations ratios, fractional contributions, etc.) 
so let us illustrate some particular cases. 

1. y = x + z where cov(x, z) = 0. This case is illustra- 
tive for population ratios, say the WC/WR ratio. The 
resulting relative error as a function of x, y is: 



u 2 x 2 y 2 xy 



(A.17) 



For the particular case of ratios of population numbers 
it is easy to show that a 2 = x = N e g(x) and thus 



1 



N eS (u) 



1 



1 



iV eff (a:) A^off(y) 



(A.18) 



y = r + z where cov(r, z) — ; but cov(r,x) ^ 0: This 
case is illustrative for ratios of line luminosities, say 
the L(WRbump)/L(H/3) ratio, where the WR has a 
contribution to L(H/3). The resulting relative error is 



^cov(r, x) 
xy 



(A.19) 



In this case, we need to know the value of cov(r, x), say 
the contribution of WR stars to the H/3 luminosity, in 
order to evaluate the error. 
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y = r + z where cov(r, z) 7^ 0, but cov(r, x) =/= 0: 
This case is illustrative for ratios of equivalent widths, 
where, for example, x is the L(H/3), r is the nebular 
continuum and z the stellar one. The resulting relative 
error is 



a x <J7 r + a z + 2 cov(r, z) 
x 2 y 2 
^cov{x, r) + cov(x, z) 
xy 



(A.20) 



which becomes Eq. A. 16. In the study of covariance 
terms, it is not necessary to worry about which arc 
the contributions if the covariance terms have been 
properly taken into account. 

A. 6. Error in a logarithm 

Last, but not least, let us assume that u = a log a; + c, 
where a and c are constants. Using Eq. |A.l| : 

_2 „2 

(A.21) 



x 2 N eS (x) 

The denormalised value of the effective number of, say, 
N c g(B), where B is the Johnson B band, gives us the ab- 
solute error in the B magnitude for the given mass trans- 
formed into stars. 

For the case of colours indices, say U-B, using the Eq. 



A.1C 
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4-b - (2-5 x log 10 (e)) : 
cov(Lu,L B ) 



Lu Lb 



V cff (Lu) V off (L B ) 



(A.22) 



Note that a N e s(U-V) value cannot be obtained due 
to the logarithm dependence in the luminosities. This sit- 
uation is general for all the colors indices. 
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